I. Web Scraping

BaseballSavant is a website dedicated to providing the public with access to historical MLB data. These data include video footage and Statcast tabular data for every pitch thrown in the MLB since 2018 and 2015, respectively. We built web scrapers to retrieve both the video source URL and pitch type for every pitch of interest. Video source URLs were used as inputs for our feature engineering pipeline.

II. Feature Engineering

A highlight in our work is the feature engineering process, termed the Video2Data pipeline. The pipeline works as follows. First, a video is downloaded from the source URL and is converted to a series of images (or frames). Second, an object detection model is used to determine the location of the pitcher in each frame. The coordinates reported by the model are subsequently used to blur the background of each image. Blurring is performed using OpenCV. The object detection model used in this step is a custom Detectron2 model (Faster R-CNN) that was trained on a self-annotated data set to specifically detect pitchers. This step is necessary for scalable and reliable feature extraction since the OpenPose pose estimation software (used in the following step) detects humans non-specifically. Third, the OpenPose pose estimation software is used on each image to extract the coordinates of 25 keypoints on the pitcher’s body. Keypoint coordinates from each frame are finally merged into a similar data structure. Example outputs generated at each step of the Video2Data pipeline are shown in Figure 1.

Figure 1. Example Video2Data pipeline outputs. (1) Video to image conversion (left). (2) Pitcher detection and background blurring (middle). (3) OpenPose pose estimation (right).

III. Data Preprocessing

Videos accessible through BaseballSavant are inherently inconsistent (i.e., in terms of video duration, camera perspective, etc). Additionally, the Video2Data pipeline produced some undesired artifacts that could be harmful for machine learning applications (e.g., missingness and pose estimation errors). As such, we implemented a data cleaning process that both detects/removes unusable observations and prepares usable observations for modeling.

The first step of our data preprocessing method addresses the inconsistent video duration problem. To crop each video to a similar range, we use pose estimation data to identify the frame at which the pitcher’s knee reaches it’s \(y_{max}\), which we refer to as \(t_{y_{max}}\) (i.e., the frame at which the left knee’s \(y\)-coordinate is maximized if the pitcher is right-handed, and vice versa). We decided to use \(t_{y_{max}}\) as a reference time point since \(t_{y_{max}}\) reliably corresponds to a common event in a pitcher’s delivery (termed the peak leg lift), irrespective of who is pitching. We then use \(t_{y_{max}}\) to determine the start and end frames, termed \(t_{start}\) and \(t_{end}\), using Equations 1-3. Briefly, \(t_{start}\) and \(t_{end}\) are computed for each observation and are used to crop each observation’s OpenPose coordinate data to a common 15 time points. An example of the frame range included for our prediction task is shown in Figure 2. Observations with \(t_{start} < 0\) are considered to have an insufficient number of frames for inclusion and are removed from consideration.

\[ n_{frames} = 15 \tag{Equations 1-3}\\ t_{end} = t_{y_{max}} + 2\\ t_{start} = t_{end} - (n_{frames}-1)\\ \]

Figure 2. An example of the cropped frame range used for pitch classification.

The second step of our data preprocessing method aims to further remove erroneous observations. Three selection criteria are used for evaluation, outside of the criteria established previously. First, we require the angle between each knee at \(t_{y_{max}}\) (the peak leg lift) to be between 40° and 90°. This selection criteria allows us to identify video footage that is taken from an atypical camera perspective, given that the standard camera perspective faces from center field to home plate. For example, in no circumstance would the angle between knees at \(t_{y_{max}}\) be between 40° and 90° if the video camera is facing from third to first base. Second, we require observations to have at least 80% of \(x\), \(y\) coordinate pairs available. Third, we perform a manual inspection using data visualizations to identify erroneous observations that were not filtered using the aforementioned selection criteria.

The third and final step of our data preprocessing method involves preparing selected observations for modeling. First, we handle missing coordinate data using the pandas time series aware interpolate() function. Second, we select 10 keypoints (or body parts) to include in downstream applications. These keypoints include the neck, the left and right hip, the left and right knee, the left and right ankle, and the wrist, elbow, and shoulder of the pitcher’s throwing arm (e.g., the right wrist, elbow, and shoulder for right-handed pitchers). Third, we normalize all \(x\) and \(y\) coordinates for a given observation to be between 0 and 1. The \(x\) coordinates for a given observation are normalized using \(x_{min}\) and \(x_{max}\) from the full set of \(x\) coordinates for that particular observation (i.e., the \(x\) values from all keypoints). The \(y\) coordinates for a given observation are similarly normalized. The aim here is to address between-video variation in pitcher location and scale. This normalization step is not to be confused with the feature scaling step that is used during modeling. Fourth, and finally, we label our observations with the pitch type outcome. Many well-cited articles dichotomize the pitch prediction task into a binary classification problem of predicting fastball versus non-fastball. In this report, we dichotomize the pitch prediction task into three classes based on Statcast labeling: fastball (FB), fastball with movement (FBwM), and off-speed (OFF). This dichotomization is presented in Table 1. Note that Tyler Glasnow’s pitch type set includes only fastballs and off-speed pitches. Thus, Tyler Glasnow’s pitch prediction task is dichotomized into a binary classification problem. The final feature set was three-dimensional per observation and took the shape of (15, 10, 2) (i.e., 15 time points, 10 body parts, and two coordinates (\(x\), \(y\))).

Model Specifications

We designed four candidate convolutional neural network (CNN) architectures for our prediction task. Architecture specifics are described in the section below. Constant model hyperparameters and training methods are specified here. All models utilize the Adam optimizer with a starting learning rate of \(1 \times 10^{-4}\), and a batch size of 16. Learning rates were reduced by a factor of \(0.5 \times 10^{-2}\) during training if the validation loss did not improve for 10 epochs. Similar dropout layers and regularization penalties were used for each model. We used a dropout rate of 0.3 in dropout layers and applied L1 and L2 regularization penalties (\(1 \times 10^{-2}\) and \(1 \times 10^{-4}\), respectively) on the kernel of each model’s last three (non-output) dense layers. The maximum number of epochs was set to 200. Training was stopped early, however, if validation loss did not decrease over the course of 50 epochs. Binary cross-entropy and categorical cross-entropy were used as the loss functions for binary and multiclass classification models, respectively. Loss functions were weighted according to outcome class weights to account for class imbalance. Binary and multiclass classification models use sigmoid and softmax activation functions in the output layer, respectively. All models were implemented using Keras (version 2.12.0) in Python (version 3.10.11).

Model Evaluation

Prior to modeling, we split the data set into a cross-validation (CV) set (80% of the data) and holdout test set (20% of the data). CV and test set sample sizes are shown for each pitcher in Table 2. Data set splits were stratified by the outcome to ensure balanced class representation across each set. We then used stratified \(k\)-fold CV with five folds to evaluate each candidate architecture. During 5-fold CV, the CV data set is split into 5 folds. Five separate models (per candidate model) are then trained using \(k-1\) (four) folds as the training data set and one fold as the validation data set. The original holdout test set is used to evaluate each of the five models. Features were standardized by removing the mean and scaling to unit variance. The mean and variance values used for standardization were learned from the training set and applied to the validation set and holdout test set during each round of 5-fold CV.

Prediction accuracy, area under the receiver operating characteristic curve (ROC-AUC), area under the precision-recall curve (AUPRC), Matthew’s correlation coefficient (MCC), and weighted F1-score were used to evaluate the predictive performance of each model. We report the mean of each performance metric (from the aforementioned 5-fold CV, using the original holdout test set for evaluation) for each model architecture. Additionally, we use one-vs-rest (OvR) ROC-AUC and AUPRC for multiclass classification evaluation since ROC-AUC and AUPRC are traditionally used for binary classification. The one-vs-rest approach binarizes the predictions into reference and non-reference classes to compute metrics (in this case ROC-AUC and AUPRC) for each class. MCC serves as a useful all-purpose metric for between-pitcher comparisons since it can be used to evaluate both binary and multiclass prediction performance.

I. CNN-2D

The CNN-2D model was originally designed as a baseline model that the more complex model architectures could be compared to. We term this model architecture CNN-2D since we reshape the 3-dimensional (3D) input shape of (15, 10, 2) to a 2-dimensional (2D) shape (10, 30) (i.e., the time steps and \(x\), \(y\) coordinates are collapsed). We additionally tried using a 2D input shape of (15, 20), whereby the body part dimension and \(x\), \(y\) coordinates are collapsed, but ultimately chose the (10, 30) configuration due to improved performance. The CNN-2D architecture consists of four convolutional layers that are connected to a series of six dense layers. Convolutional layers consist of 2D convolutions, using kernels of size (3, 1), followed by batch normalization and rectified linear unit (ReLU) activation functions. We found that (3, 1) kernels outperformed (3, 3) kernels in terms of ROC-AUC and accuracy. An increasing number of filters (8, 16, 32, and 64 filters) were used for consecutive convolutional layers. The output of the convolutional block is connected to a series of six fully connected dense layers. The first five dense layers incorporate 256 (with dropout), 128 (with dropout), 64 (with regularization), 32 (with regularization), and 16 (with regularization) nodes and ReLU activations. The final dense layer (the output layer) consists of one node when the sigmoid activation function is used and three nodes when the softmax activation function is used.

II. CNN-3D

The CNN-3D model architecture strongly resembles the CNN-2D model architecture, although the 3D input shape is preserved (hence CNN-3D). Distinct from the CNN-2D model architecture, convolutional layers use 3D convolutions with kernels of size (3, 3, 1). We were restricted to using this kernel size due to the size of the 3D input. All other model architecture details are consistent with the CNN-2D model architecture.

III. CNN-3D-LSTM

The CNN-3D-LSTM model architecture builds off of the CNN-3D model architecture by incorporating a 256-unit long short-term memory (LSTM) layer between the convolutional block and dense layers. All other model architecture details are identical to the CNN-3D model architecture.

IV. Branched CNN

The Branched CNN model architecture was inspired from an article that uses a convolutional neural network to predict mortality from 12-lead electrocardiogram voltage data⁸. Briefly, the Branched CNN model architecture consists of 10 separate convolutional branches that independently extract features from each keypoint (or body part). The features extracted by each branch are concatenated and fed to a series of fully connected dense layers that are identical to all other model architectures. Similar to the convolutional block used in the CNN-2D architecture, each branch consists of four convolutional layers that use 2D convolutions and kernels of size (3, 1). However, ReLU activations are (unconventionally) used directly after convolutions (and before batch normalization). We found that that the Conv2D-BatchNormalization-ReLU configuration used in other model architectures led to inconsistent and diminished model performance. Additionally, global average pooling is independently applied to the output of each convolutional branch prior to concatenation.